#include <stdio.h>
#include <stdlib.h>

#define ABORT(msg) { fprintf(stderr, "%s\n", msg); exit(1); }

typedef struct {
    int vertex;
    int degree;
} QUEUE;
int cmp(const void *a, const void *b) { return ((QUEUE *)a)->degree - ((QUEUE *)b)->degree; }

typedef struct {
    int removed;    // 从图中移除的节点
    int visited;    // 被访问
    int degree;     // 节点度数
    int *neighbors; // 邻接链表
} VERTEX;

typedef struct {
    VERTEX *vertexs;

    // buffer
    int *adjency;
    int *depth;
} GRAPH;

GRAPH* construct_graph(const int n, const int *row, const int *col, const int threshold) {
    VERTEX *vertexs = (VERTEX *)calloc(n, sizeof(VERTEX));

    int nnz = col[n];

    int x, y;

    // 度数统计
    for (int j = 0; j < n; j++) {
        for (int i = col[j]; i < col[j + 1]; i++) {
            x = row[i], y = j;
            vertexs[x].degree++;
            vertexs[y].degree++;
        }
    }

    // 标记关键节点
    for (int i = 0; i < n; i++) {
        vertexs[i].removed = vertexs[i].degree > threshold;
        vertexs[i].degree = 0;
    }

    // 剪枝，去掉关键节点
    for (int j = 0; j < n; j++) {
        for (int i = col[j]; i < col[j + 1]; i++) {
            x = row[i], y = j;
            if (!vertexs[x].removed && !vertexs[y].removed) {
                vertexs[x].degree++;
                vertexs[y].degree++;
            }
        }
    }

    // 分配邻接链表存储空间
    int size = 0;
    for (int i = 0; i < n; i++) { size += vertexs[i].degree; }
 
    int *adjency = (int *)malloc(size * sizeof(int));

    int offset = 0;
    for (int i = 0; i < n; i++) {
        vertexs[i].neighbors = &adjency[offset];
        offset += vertexs[i].degree;
    }
    
    // 构建邻接链表
    for (int j = 0; j < n; j++) {
        for (int i = col[j]; i < col[j + 1]; i++) {
            x = row[i], y = j;
            if (!vertexs[x].removed && !vertexs[y].removed) {
                vertexs[x].neighbors[ vertexs[x].visited++ ] = y;
                vertexs[y].neighbors[ vertexs[y].visited++ ] = x;
            }
        }
    }
    for (int i = 0; i < n; i++) vertexs[i].visited = 0;   

    GRAPH *graph = (GRAPH *)malloc(sizeof(GRAPH));
    graph->vertexs = vertexs;
    graph->adjency = adjency;
    graph->depth = (int *)malloc(n * sizeof(int));
    return graph;
}

void destroy_graph(GRAPH *graph) {
    if (graph) {
        free(graph->vertexs);
        free(graph->adjency);
        free(graph->depth);

        free(graph);
    }
}

int neighbor_to_queue(const int n, GRAPH *graph, const int vertex, QUEUE *queue, int offset) {
    int len = 0;
    VERTEX *vertexs = graph->vertexs;

    for (int i = 0; i < vertexs[vertex].degree; i++) {
        int x = vertexs[vertex].neighbors[i];
        if (!vertexs[x].visited && !vertexs[x].removed) {
            vertexs[x].visited = 1;
            queue[offset + len].vertex = x;
            queue[offset + len].degree = vertexs[x].degree;
            len++;
        }
    }

    return len;
}

void check(const int n, const int *perm) {
    int *a = (int *)calloc(n, sizeof(int));
    for (int i = 0; i < n; i++) {
        if (perm[i] >= n)  ABORT("cuthill mckee error");
        a[ perm[i] ] = 1;
    }

    for (int i = 0; i < n; i++) {
        if (a[i] != 1) ABORT("cuthill mckee error");
    }
    printf("check OK\n");
    free(a);
}

// 返回该子图的顶点个数
int bfs2(int n, GRAPH *graph, int *start_vertex, int *deepest, int visited, QUEUE *queue) {
    VERTEX *vertexs = graph->vertexs;
    int *depth = graph->depth;

    int left = 0, right = 0;
    queue[right].vertex = *start_vertex;
    queue[right].degree = vertexs[*start_vertex].degree;
    vertexs[*start_vertex].visited = 1;
    depth[right] = 1;
    right++;

    while (left != right) {
        int index = queue[left].vertex;

        int size = neighbor_to_queue(n, graph, index, queue, right);
        // qsort(&queue[right], size, sizeof(QUEUE), cmp);
        for (int i = 0; i < size; i++) {
            depth[right + i] = depth[left] + 1;
        }
        left++;
        right += size;
    }

    int index = right - 1;
    *deepest = depth[index];

    for (int i = right - 1; i >= 0 && depth[i] == *deepest; i--)
        index = queue[index].degree < queue[i].degree ? index : i;

    *start_vertex = queue[index].vertex;

    if (!visited)
        for (int i = 0; i < right; i++) vertexs[ queue[i].vertex ].visited = visited;

    return right;
}

int get_connected_compontend(int n, GRAPH *graph, int *vertexs, int *connected) {
    int len = 0;
    int deepest;
    QUEUE *queue = (QUEUE *)malloc(n * sizeof(QUEUE));
    for (int i = 0; i < n; i++) {
        if (!graph->vertexs[i].removed && !graph->vertexs[i].visited) {
            vertexs[len] = i;
            connected[len] = bfs2(n, graph, &vertexs[len], &deepest, 1, queue);
            len++;
        }
    }
    free(queue);
    for (int i = 0; i < n; i++) graph->vertexs[i].visited = 0;
    return len;
}

/*
目的使得遍历数越深越好，每层的节点数量就会更少
https://zhuanlan.zhihu.com/p/20506574
1. 任意初始a0，层序遍历，计算层数d0
2. 取最后一层度最低节点a1，层序遍历，计算层数d1
3. 若d1 > d0, a1作为初始节点，重复，否则输出a1
*/
void pseudo_peripheral2(int n, GRAPH *graph, int *vertex, QUEUE *queue) {
    int a0, d0;
    int a1 = *vertex, d1;
    bfs2(n, graph, &a1, &d1, 0, queue);

    do {
        a0 = a1, d0 = d1;
        bfs2(n, graph, &a1, &d1, 0, queue);
    } while (d1 > d0);

    *vertex = a1;
    if (d1 == d0) *vertex = graph->vertexs[a0].degree < graph->vertexs[a1].degree ? a0 : a1;
}

void update_start_vertexs(int n, GRAPH *graph, int vertexs_size, int *vertexs, int *connected) {
    QUEUE *queue = (QUEUE *)malloc(n * sizeof(QUEUE));
    for (int i = 0; i < vertexs_size; i++) {
        pseudo_peripheral2(n, graph, &vertexs[i], queue);
    }
    free(queue);
    for (int i = 0; i < n; i++) graph->vertexs[i].visited = 0;
}

#if 0
void sort_by_degree(int n, QUEUE *queue) {
    int i, j;
    QUEUE key;
    for (int i = 1; i < n; i++) {
        key.vertex = queue[i].vertex;
        key.degree = queue[i].degree;
        j = i - 1;
        while (j >= 0 && queue[j].degree > key.degree) {
            queue[j + 1].vertex = queue[j].vertex;
            queue[j + 1].degree = queue[j].degree;
            j--;
        }
        queue[j + 1].vertex = key.vertex;
        queue[j + 1].degree = key.degree;
    }
}
#endif

// TODO:如果有连通分量大于1，则剖分。
// Only support CSC
void partial_reverse_cuthill_mckee(const int n, const int *row, const int *col, int *perm, int threshold) {
    // printf("t = %d\n", threshold);

    // 构建图
    GRAPH *graph = construct_graph(n, row, col, threshold);

    // 获取关键节点序列
    QUEUE *queue = (QUEUE *)malloc(n * sizeof(QUEUE));
    int left = 0, right = 0;
    for (int i = 0; i < n; i++) {
        if (graph->vertexs[i].removed) queue[left++].vertex = i;
    }
    right = left;

    // 获取所有的连通分量
    int *vertexs = (int *)malloc(2 * n * sizeof(int));
    int *connected = vertexs + n;
    int vertexs_size = get_connected_compontend(n, graph, vertexs, connected);

    // 更新初始节点
    update_start_vertexs(n, graph, vertexs_size, vertexs, connected);

    // 遍历
    for (int i = 0; i < vertexs_size; i++) {
        int start = vertexs[i];
        queue[right++].vertex = start;
        graph->vertexs[start].visited = 1;
        graph->vertexs[start].removed = 1;

        while (left != right) {
            int index = queue[left++].vertex;

            int len = neighbor_to_queue(n, graph, index, queue, right);
            qsort(&queue[right], len, sizeof(QUEUE), cmp);
            right += len;
            for (int i = 0; i < graph->vertexs[index].degree; i++) {
                graph->vertexs[ graph->vertexs[index].neighbors[i] ].removed = 1;
            }
        }

    }
    for (int i = 0; i < n; i++) perm[i] = queue[i].vertex;
    destroy_graph(graph);

#ifdef DEBUG
    check(n, perm);
#endif

    // reverse
    int i = 0, j = n - 1;
    while (i < j) {
        int t = perm[i];
        perm[i++] = perm[j];
        perm[j--] = t;
    }

    free(vertexs);
    free(queue);
}

